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Using Wannier function-based interpolation techniques, we present compelling numerical evidence 
for the presence of a saddle-point van Hove singularity at the F point near the phosphorene Fermi 
energy. We show that in proximity of the van Hove singularity the spin susceptibility presents the 
logarithmic temperature dependence typical of Liftshitz phase transitions. Furthermore, we demon¬ 
strate that the critical temperature for the ferromagnetic transition can be greatly increased (up to 
0.05 K) if strain along the zigzag ridges is applied. Although the ferromagnetic state would be very 
difficult to experimentally reach, the logarithmic temperature behaviour of the spin susceptibility 
due to the van Hove singularity is found to persist at much higher temperatures (up to ~97 K). 
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I. INTRODUCTION 

Saddle-point van Hove singularities (VHSs) originate 
from saddle points in the band structure, around which 
the band curvature has opposite signs along two orthog¬ 
onal directions. In two dimensions, the density of states 
(DOS) diverges at the VHS, and therefore arbitrary weak 
interactions can produce large effects in the electronic be¬ 
haviour, giving rise to instabilities in many aspects such 
as charge, spin, and/or pairing susceptibilities. Once 
the Fermi energy ajmroaches a VHS, ferromagnetism,!^ 
antiferromagnetism,® and/or superconductivitj®® can 
be substantially enhanced. 

The VHS is a topological critical point of the Fermi 
surface, across which the quantum Lifshitz phase tran¬ 
sition takes place.®tni The Lifshitz transition for non¬ 
interacting systems is continuous and does not break 
symmetry. For interacting systems, however, the Lif¬ 
shitz transition may beco me discontinuous and accom¬ 
pany symmetry breaking.®!^ In cuprates, Hall coeffi¬ 
cient measurements provid e evidence for the Fermi sur¬ 
face topology change .!^2114| xhe Lifshitz transition is also 
proposed to change the Fermi liquid into the marginal 
Fermi liquid,!^^ and the VHS is thus argued to be respon¬ 
sible for the linear temperature (T) dependence of resis¬ 
tivity and the T-independent thermopowei!^^ observed in 
this regimeMoreover, in the so-called “van Hove 
scenario”, the presence of a VHS near the Fermi energy 
is argued to play a major role in the high-Tc supercon¬ 
ductivity of cuprates! 23|2i | Given the strong influence of 
VHSs on the properties of materials, it is important to 
identify the presence and understand the role of these 
singularities, especially for technologically promising low¬ 
dimensional m ateria ls like phosphorene. 

Phosphorene,!22123l ^ single layer of black phosphorus, 
is the most recent addition to the growing family of two 
dimensional (2D) materials. It is a semiconductor with 
high potential for applications in electronic and optoelec¬ 
tronic devices.!^ Despite the relative infancy of the held, 
few-layer phosphorene held ehect transistors exhibit very 


high on-off current ratio^^HHl (exceeding 10^) and am- 
bipolar behaviour,!^ together with the highest hole mo¬ 
bility eve^dOOO cm^/Vs) for a 2D material apart from 
graphene.!23 Phosphorene’s pliable waved structure also 


allows for strain engineering of both ehective masses and 
bandgaps.!^ Strain can even induce a semiconductor to 
metal transition.!^ 

In this article, we show that a VHS is present at the 
phosphorene Fermi energy, and we investigate the con¬ 
sequent ferromagnetic instability in both the unstrained 
and strained cases. 

The article is organized as follows: after introducing 
the computational methodology, in Sec. |HI A| we present 
the electronic structure of phosphorene near the VHS, in 
Sec. HIB we study the ferromagnetic instability (with¬ 
out strain), and finally in Sec. HIC we investigate the 
effect of strain on the critical temperature of the fer¬ 
romagnetic transition. 


II. COMPUTATIONAL DETAILS 

The calculations involve the following three steps. 

(i) A density functional theory (DFT) calculation is per¬ 
formed with a plane wave basis set, as implemented in the 
Quantum ESPRESSO package.!® We use the PBEsol 
functionapi! for the exchange and correlation energy. A 
plane wave basis set with a kinetic energy cutoff of 70 Ry 
(280 Ry) is used to represent the electronic wave func¬ 
tion (charge density). The core electrons are described 
via the projected-augmented wave (PAW)(® method; 12.9 
A of vacuum are added in the direction normal to the 
monolayer to avoid spurious interactions between peri¬ 
odic replicas. Both lattice parameters and atomic posi¬ 
tions are relaxed until the forces on each atom are less 
than 10“3 eV/A and the pressure is less than 1 kbar. 
After lattice relaxation, the phosphorene crystal param¬ 
eters are a2,=3.28 A and ay=4.44 A, in agreement with 
a previous study.!® xhe optimized configuration of the 
phosphorene monolayer is presented in Fig. [0a). In this 
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DFT calculation, the Brillouin Zone (BZ) is sampled us¬ 
ing a F-centered 60x48x1 Monkhorst-Pack (MP) gridP^ 
This calculation will serve as a benchmark for the Wan- 
nier interpolation of the band structure [Fig. and 

Fig- 3c)]- 

(ii) Using the self-consistent charge density obtained 
from step (i), we evaluate the required input quantities 
for the Wannier calculation (energy eigenvalues, over¬ 
lap matrices and projections^S) on a relatively coarse 
10x8x1 mesh for the unstrained case (Secs. Ill A and 


IIIB), and a 30x8x1 mesh for the strained case (Sec. 


IIIC). These fc-point meshes are fine enough to pro¬ 
vide converged Wannier functions. The calculations are 
performed with Quantum ESPRESSO and its post¬ 
processing subroutine pw2wannier90. 

The aim of the two previous step s is t o obtain the 
maximally localized Wannier function^^HSH (MLWFs) to 
be later used for the very dense fc-point sampling around 
the VHS. 


(iii) With the energy eigenvalues, overlap matrices and 
projections obtained from step (ii), we construct the ML¬ 
WFs according to the procedure presented in Ref. [5^ and 
Ref. |36l The resulting Wannier functions consist of three 
p-orbitals centered on each P atom, leading to the wan- 
nierization of 6 valence and 6 conduction bands (there 
are four P atoms in the phosphorene unit cell). 

One of the main advantages of the maximally localized 
Wannier representation of the DFT orbitals is that quan¬ 
tities calculated on a coarse reciprocal-space grid can be 
used to interpolate on a much finer grid with low compu¬ 
tational cost. The Wannier interpolation is particularly 
useful when a fine BZ sampling is required to converge 
the quantity of interest. In this work, such quantities are 
the DOS and the valence band in a small region around 
the VHS. For the DOS calculation, an extremely dense 
Wannier interpolated mesh of 4-10"^x3.2-10^x 1 - corre¬ 
sponding to ^1.3 billion /c-points in the BZ - is used to 
capture the sharp peak in the DOS due to the VHS. 
We use a smearing of 7x10“^ eV. Similarly, a Wan¬ 
nier interpolation with a reciprocal lattice spacing of 
Akx^y=‘ 2 - 10 ~^ X corresponding to a 500x500x1 

MP grid is employed for the contour plot of the valence 
band around the VHS. All MLWF calculations are per¬ 
formed with the Wannier90 package.!^ 

To further validate the PBEsol results, we have 
performed additional calculations with the local LDA 
functional,!^ the semilocal PBE functionaP^, the 
screene d-hyb rid HSE06 functional,!^ and the GW 
methocpIBll order to elucidate the position of the va¬ 
lence band maximum. In particular, the HSE06 func¬ 
tional and the GW method are known to provide a more 
accurate description of the electronic properties of semi¬ 
conductors and i nsulators than local or semilocal density 

functionals.EHBll 


In LDA and PBE calculations, atomic positions and 
lattice parameters are relaxed till the forces are less than 
10“^ eV/A and the pressure less than 1 kbar. The BZ is 
sampled with a F-centered 60x48x1 MP grid as in the 


case of PBEsol calculations. The HSE06 corrections are 
instead calculated self-consistently using the PBE relaxed 
lattice parameters and atomic positions, together with a 
12x12x1 MP grid. The band structure is then obtained 
using the derived Wannier functions in a similar fashion 
to the PBEsol calculations outlined above. 

The GW calculations are performed in two steps. 
First, atomic positions and relaxed lattice geometries are 
calculated with the PBE functional and norm-conserving 
Troullier-Martins pseudopotentials.!^ Then the GW cor¬ 
rections are computed following the method proposed by 
Hybertsen and Louie.!^ We include 138 bands in the eval¬ 
uation of the dielectric matrix and the self-energy, with a 
cutoff of 4 Ry for the dielectric matrix; convergence was 
checked including up to 384 bands. A supercell of 20 A in 
the direction perpendicular to the monolayer and a slab- 
truncation potentiaP^ are used in these GW calculations 
to avoid spurious interactions with periodic replicas of 
the system. A fine MP grid is employed in the direction 
of the VHS (100x8x1) to distinguish the top of the va¬ 
lence band from the F point, thus providing evidence for 
the presence of the VHS. All GW calculations are Mr- 
formed with the plane-wave based ABINIT package.!^ 


III. RESULTS AND DISCUSSION 


A. Electronic band structure and van Hove 
singularity 


(a) (b) 



Figure 1: (Color online) (a) Crystal structure of 
phosphorene and its projections on the y-z, x-y and x-z 
planes (b) DFT-PBEsol electronic band structure (solid 
black line) and its Wannier interpolation (dashed red line). 
The Brillouin zone is also shown (c) Detail of the electronic 
band structure and (d) DOS in a small region near the VHS. 

The energy at the VHS is set to zero. 

Phosphorene is a semiconductor with a relatively large 
bandgap that is underestimated (0.72 eV) at the PBEsol 
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level (a well-known deficiency of local and semilocal DFT 
functional s^, and enlarged at 1.6 — 2.0 eV when GW 
correction's^ are included .SStSS The electronic band 
structure of phosphorene, calculated with the PBEsol 
functional, is shown in Fig. Bb). The black solid line 
represents a standard plane wave DFT calculation, while 
the red dotted line is the band structure obtained through 
a Wannier interpolation. Our Wannier interpolation is 
very accurate, over a broad range of energies. In particu¬ 
lar, at F point, the DFT band structure and the Wannier 
interpolation differ by less than 10“^ eV. 

In agreement with recent studies,l2SIIIIIl find that 
the top of the valence band is slightly away from the F 
point for the LDA, PBE and PBEsol functionals. Using 
these three DFT functionals, we consistently find that the 
top of the valence band is displaced from F along the F-X 
direction, which is the direction along the phosphorene 
zigzag ridges [see Fig. Qa)]. The detailed symmetry 
analysis presented in Ref. [52] attributes the absence of 
direct bandgap to the counteracting effects (in the k • p 
approximatiorP^ of states of different symmetries on the 
valence band around the zone center. 

To further validate these results, we have carried out 
calculations using the screened hybrid HSE06 functional 
and the GW approximation, which are known to improve 
the description provided by local or semilocal DFT func¬ 
tionals (such as LDA, PBE and PBEsol), not only regard¬ 
ing bandgaps, but also concer ning the band dispersion in 
semiconductors and insulators .HSElISlHSll 

To quantitatively characterize the valence band maxi¬ 
mum, we define kniax=(fcmaxi 0) as the wavector at which 
the valence band has a maximum, and Emax as the differ¬ 
ence in energy between the valence band maximum and 
the value of the valence band at the P point. The re¬ 
sults obtained with various computational methods and 
different strains are reported in Table]^ As mentioned be¬ 
fore, in absence of strain, LDA, PBE and PBEsol gives 
a valence band top slightly away from the P point, with 
Aniax ranging approximately from 1 to 12 meV. With 
the HSE06 functional, the valence band top is slightly 
misplaced from F; however, the calculated value of Emax 
is so small (0.05 meV) that it can be considered zero. 
Also the GW method - in the absence of strain - predict 
phosphorene to be a direct bandgap semiconductor. 

We then apply strain along the cc-direction, changing 
the lattice parameter Ux to be ax{l -l-a^str), where a pos¬ 
itive (negative) value of Xstr indicates tensile (compres¬ 
sive) strain. The results are shown in Table[lj Application 
of compressive strain moves the valence band top away 
from the F point in all computational methods. In par¬ 
ticular, with the HSE06 functional and the GW method, 
the top of the valence band is displaced from the F point 
with a 2% strain; larger strains monotonically increase 
both fcniax and ifmax with ifmax~ 54 meV for a strain of 
8% according to the GW method. In contrast, tensile 
strain moves the top of the valence band towards the F 
point, and eventually removes the VHS singularity. 

Having established that a VHS near the phosphorene 


Fermi energy is either present or can be strain-induced 
using a wide range of electronic structure descriptions, 
hereafter we consider as an explanatory example the case 
of the PBEsol functional. Other functionals and the GW 
method are expected to yield similar general results. 

A magnification of the valence band maximum is 
shown in Fig. [^c). From Fig. [^c), we notice that the 
valence band has a saddle point at F. In the reciprocal 
space neighbourhood of this point, the principal curva¬ 
ture is electron-like along the F-X path [from the left in 
Fig. 0 c)], while it is hole-like in the F-Y path [from 
the right in Fig. Qc)]. Thus, at the F point there is a 
crossover from electron-like to hole-like conduction that 
originates at the VHS. The DOS, calculated on an ultra- 
fine grid of ^1.3 billion fc-points, is shown in Fig. [^d). It 
exhibits a divergent behaviour at the energy position of 
the VHS, as expected for a 2D lattice. In contrast to the 
saddle point behaviour at F, the valence band has a max¬ 
imum at kmax and therefore the DOS shows a step-like 
drop to zero at this point 


(a) 



(b) 



Figure 2: (Color online) (a) 3D plot of the phosphorene 
valence band around the VHS (F point) (PBEsol 
functional). Cut-x(y) indicates the r-X(Y) path used in the 
band structure calculation, (b) 2D contour plot of the 
valence band in a larger region around F. The contour lines 
are drawn at 2 meV intervals. The energy at the VHS is set 
to zero in both plots. 

Three dimensional (3D) and 2D plots of the phospho¬ 
rene valence band in the neighbourhood of the VHS (F 
point) are depicted in Fig. [^a) and Bb), respectively. 
The electron-like dispersion along the F-X path [Cut-x 
in Fig. Ba)] and the hole-like dispersion in the F-Y path 
[Cut-y in Fig. B^)] evident. The VHS has indeed 
the topology of a 3D-saddle point. Moreover, the valence 
band is anisotropic at the F point, with strong disper- 
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11.6 
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0.067 
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0.000 

0.0 
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0.000 
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Table I: Position of the valence band maximum, fcmaxj and its energy, -Emaxj relative to the P point [see FigQc)] 
calculated with different computational methods and various strains, ccgtr, along the x-direction of the phosphorene 
lattice. When kjnax=0 (and thus i?max=0) the top of the valence band coincides with the P point, and therefore the 

VHS is not present. 


sion along P-Y [armchair direction, see Pig. [^a)], while 
it is nearly flat on P-X [zigzag direction, see Pig. [^a)]. 
The large difference in magnitude of the effective masses 
along the two directions gives to the VHS an extended 
structure, as shown in Fig. [^b). From a fitting of the lo¬ 
cal curvature of the valence band around the P point, we 
obtain maj/ruy ~27, where rux and ruy are the effective 
masses on the P-X and P-Y path, respectively. 

In the limit of infinite mass in one direction (i.e. flat 
band in one direction, vanishing curvature) the sad¬ 
dle point becomes extended, giving rise to a so-called 
extended VHS (EVHS). EVHSs have been experimen¬ 
tally observed in doped graph enj^ and in some layered 
cuprate superconductorsjn 2D materials, the DOS 
is known to diverge logarithmically at the VHS, while 
in an EVHS the energy dispersion is quasi-one dimen¬ 
sional, and the DOS has a much stronger square-root 
divergence.!^ Therefore, due to the anisotropy of the 
phosphorene band structure, the VHS has an extended 
character that might amplify its effects on the material 
properties. 


B. Ferromagnetic instability 


As mentioned in the introduction, the presence of a 
VHS at the Fermi energy can create ferromagnetic, anti¬ 
ferromagnetic or superconducting instabilities. In con¬ 
trast to cuprates where the VHS points are at (tt, 0) 
and (0, tt), in phosphorene the VHS point is at P and 
therefore we can exclude antiferromagnetism since no 
inter-VHS scattering can induce this instability. Fur¬ 
thermore, for highly anisotropic masses (see Sec. HI A I, 
"nix/my 3> 1, similar to the t — t' Hubbard model with 
large t'/t (> 0.276), th e ferro magnetic instability will win 
over other instabilities! ^ * ^!^"^ ! As a result, we can omit also 
superconductivity and consider only ferromagnetism. 

The extremely fine structure of the VHS in phospho¬ 
rene requires a very high resolution calculation of the 
band structure. To make the calculation accessible, in¬ 


stead of using the band structure from the Wannier in¬ 
terpolation, we approximate it here by an analytic single¬ 
band model. Consistent with Fig. Qc) and Fig. the 
low-energy physics in the neighbourhood of the VHS can 
be described by 

E{kx,ky) = \akl-\pkt-\a'kl ( 1 ) 


which characterizes the saddle point at P (opposite signed 
band masses along kx and ky) and band inflection along 
kx- As in the previous section, the VHS energy at P, 
EvhS) is set to zero while the band maximum at fcmax is 
Eniax- To fit the DFT-PBEsol calculations, the band 
parameters follow the relations: a'/a=mx/my=27.02, 
A/a//3=|/cniax|=0.104 and a^/4/3=Einax=6.6 meV. 

This simple model captures the energy dispersion be¬ 
haviour near the VHS and, since the parameters of the 
model are determined directly from the DPT calculation, 
it allows us to investigate the magnetic instability quan¬ 
titatively. 

Obviously, a ferromagnetic instability can take place 
only in metallic or semimetallic systems, and therefore 
some amount of doping is required for phosphorene to 
exhibit metallic behaviour. To study the effects arising 
from the presence of the VHS, hereafter we thus assume 
the Fermi energy to be exactly at the VHS, Eyns = 0, 
unless otherwise stated. In the case of the PBEsol func¬ 
tional, the amount of (hole) doping necessary to reach the 
VHS (from the top of the valence band E^^ax) is found to 
be approximately 4.2x10“^ electrons per unit cell, cor¬ 
responding to a doping concentration of 1.4x10^^ cm“^ 
for each spin. 

Given the energy dispersion in Eq. Q, it is possible to 
derive an exact analytical expression for the DOS, N{E) 
(see Appendix [A| for a complete derivation): 


N{E) = 


ot' 

2 CLxCL 


K 


l-p- 


/3a' 7r2|fc_| 
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i-pph 


if A > 0 
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( 2 ) 
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where we have defined 


k+ = fcn 




P^ = 


kt 


(3) 


and K(k) is the complete elliptic integral of the first kind. 
The quantities and ay are the phosphorene lattice pa¬ 
rameters, as defined in Sec. [ITj 

Since we are mainly interested in the behaviour of the 
DOS at the VHS, we take the limit £1—in Eq. ([^ to 
obtain (see Appendix [A|) 

N{E ^ 0±) = [ln(£;^ax/A) + 0(1)] (4) 

ZTr^Vao! 

The model of Eq. Q therefore exhibits a logarithmi¬ 
cally divergent DOS and it is proportional to the geo¬ 
metric mean mass, y/nEjfhy oc l/^/c^Q;^ so in the limit of 
small energies we can approximate the DOS as 

N{E)k Noln{A/E) (5) 


where No= ^ °‘^^^ =0.0588 eV and A is an energy 
cutoff of the order of Amax- In this approximation, the 
DOS only includes contributions from states around T. 
This is fully justified since the behaviour of the DOS 
near the VHS is obviously governed by its divergence 
at E=0 (r point), and therefore finite (not diverging) 
contributions from other regions in the Brillouin zone 
can be neglected. 

With the logarithmic DOS one can derive (see Ap¬ 
pendix 1^ an expression for the bare spin susceptibility 
that shows a dependence on the logarithm of the inverse 
temperature, 


f Sk dn ,( Et ) 


( 6 ) 

(7) 


where np is the Eermi distribution, Ek are the energy 
(Kohn-Sham) eigenvalues, and wd is a fitting constant of 
the order of Ej^six- We also set the Boltzmann constant 
kp to unity. The logarithmic divergence at low tempera¬ 
tures [Eq. ([7)] is confirmed by explicit calculation of the 
spin susceptibility using Eq. (§, as shown in Fig. |^a). 

Notably, the logarithmic behavior is present only when 
the temperature is lower than T *, which is defined as the 
temperature above which the spin susceptibility starts 
deviating from the logarithmic behaviour. Thus, above 
T* ^17 K (corresponding to an energy scale of 1.5 meV), 
we observe deviation of the susceptibility from the log¬ 
arithmic law as seen in the insert of Fig. |^a). This 
temperature T* is related to the energy scale Emax- We 
have checked this apparent relationship between T* and 
Aniax by comparing susceptibilities for different band pa¬ 
rameters /3 (and thus E^i^x) in Eq. Q, and we indeed 
see proportionality between T* and E^ax (not shown). It 
is also found that the susceptibility increases with E^^x 
as expected from the energy cutoff and fitting constant 
dependence in Eqs.([^ and Q respectively. 


(a) 




Figure 3: (Color online) (a) Temperature dependence of the 
bare spin susceptibility y calculated directly from Eq. ([^ 
(blue dots) or approximated with the logarithmic divergence 
in Eq. 0 (dashed red line). The low-temperature behaviour 
for T<r*~1.5 meV is seen in the inset to follow the 
logarithmic law. The dashed red line is fitted to Eq. 0 
with ujn =0.5004 eV. The Fermi level is set to Eyns^O. (b) 
Spin susceptibilities y different Fermi energies /r. Away 
from the VHS, the low-temperature logarithmic behaviour 
stops at T ~ and turns into Pauli susceptibility. 


Next, we examine the effect of doping on the ferromag¬ 
netic instability, considering various chemical potential 
shifts fi. The susceptibility was calculated numerically 
and the results are presented in Fig. [^b), in which the 
Fermi energy is shifted to ^ above Ayns (the energy of 
the VHS, or T). Fig. [lb) shows that even away from 
the VHS point, the logarithmic-T behaviour of the sus¬ 
ceptibility is still preserved for T<T*. However, for each 
value of fjL, we see that the logarithmic increase of y with 
decreasing T stops at /r, below which the susceptibility 
become constant suggesting a transition to Pauli param¬ 
agnetism at low temperatures. This different behaviour 
for large and small T (with respect to /r) can thus be 
understood directly from the expression of the bare spin 
susceptibility as outlined in Appendix B. 

It is in fact possible to obtain analytical estimates for 
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X in both regimes (please refer to Appendix for a com¬ 
plete derivation). For the susceptibility has the 

form 


where the momenta ki, k 2 and ka are in the neighbour¬ 
hood of the Fermi surface and the interaction strength 


X (T » /r) « TVo In (ujd/T) cosh'^ (^/2T). (8) 

We observe the logarithmic-T behaviour, typical of 
Liftshitz phase transitions. Moreover, for ^/T—>-0, 
cosh“^(/x/2T)—and therefore, in this limit, x for the 
doped system has precisely the same behaviour as in the 
undoped case, as confirmed by the numerical results pre¬ 
sented in Fig. I^b). In contrast, for the suscepti¬ 

bility is found to be independent of T: 

X (T < « A^o ln(A/^) (9) 

where A is an energy cutoff A<i?i„ax- This saturation 
of X agrees well with the numerical results presented in 
Fig. ib), and it originates from the infrared cutoff of 
the excitations due to the shifted thermal distribution 
(see Appendix . 

Now we estimate the ferromagnetic transition tem¬ 
perature. Let us assume a Hubbard interaction of 
strength U between intra-orbital spins. According to the 
Stoner criterion,!^ the magnetic transition occurs when 
C/i,x(T)=l. Here 11,^, which is defined by U times the av¬ 
erage weight Wi, at the Fermi energy for a particular or¬ 
bital z/, is regarded as the effective interaction of orbital 
u. From the Stoner criterion, the critical temperature 
follows the BCS forrn,!^ 


Tc = ud exp{-l/NoU^) (10) 

where the geometric mean mass, appearing in Nq as out¬ 
lined above, determines the DOS at the Fermi energy. 
With this BCS formula, one can obtain the magnetiza¬ 
tion directly. 

The effective interaction 14ff can be evaluated using 
the Kohn-Sham orbitals from the DFT calculation. Let 
us define orbital operators z/'m and band operators 
The relation between them is a unitary transformation 
^m(k)=X)i/ Am,i/(k)0^(k) where A(k) is the unitary ma¬ 
trix that diagonalizes the Bloch Hamiltonian. The Hub¬ 
bard onsite (intra-orbital) interaction is 

Hu = uY, (n) 

R,?71 

where R is the real space lattice vector and V'mtd are the 
Kohn-Sham spin-orbitals. 

Since only the valence band (VB) is included in our 
low-energy model [Eq. 01. we include only the intra¬ 
band scattering terms from the Hubbard model. More¬ 
over, at T=0, only states from the Fermi surface con¬ 
tribute to the susceptibility. After these considerations, 
the effective interaction Veff for z^=VB is 

Kff Y. + k3 - k2), 

ki.k2,k3 

(12) 


[/, = C7W, = C7/^|A^,,(k)|M (13) 

\ m / FS 

is averaged on the Fermi surface in the sense that the 
momentum dependence can be neglected since the Fermi 
surface around the VHS is small. 

From the DFT results, we obtain an average weight, 
HA, for the contributing orbitals of about 0.2. This or¬ 
bital weight significantly reduces the critical tempera¬ 
ture. For example, using the criterion C/^x(T'c)=0.8 and 
at C/=4 eV, the critical temperature Tc for ferromag¬ 
netism is only about 4 /rK. 

Doping can destroy ferromagnetism even at zero 
temperature when AAC/j^ ln(a;£)//i)<l. Although inter¬ 
orbital interactions might slightly enhance Tc, the Stoner 
criterion applied to the bare susceptibility typically over¬ 
estimates the critical temperature since particle-particle 
correlations would give large corrections to the self- 
energy.lS^lt^ As a result, this ferromagnetic state would 
be difficult to reach. 


C. Effect of strain on the van Hove singularity and 
on the critical temperature 


Strain can have a large effect on phosphorene’s pliable 
waved structure, and therefore it represents a natural 
way to tune the band parameters of the VHS, in order to 
increase Tc- In particular, from Eq. (10) we notice that. 


for a fixed effective interaction, the critical temperature 
can be varied in two ways. One way is to increase lod 
(or equivalently Amax, see Sec. HIB), which will cause 


a linear increase in Tc- The second and more prominent 
way is to increase AA, which will result in an exponential 
increase in T^. This can be accomplished, for instance, by 
reducing the dispersion in both x and y directions near 
the F point [see Eqs. ([^ and ([^]. 

We find that strain along the armchair direction [y- 
axis in Fig. [^a)] does not alter significantly the critical 
temperature, with a modest ninefold increase in Tc (~36 
/iK) for a tensile strain of 3%. 

The situation, however, is significantly different for 
strain along the zigzag ridge direction [x-axis in Fig. 
I^a)]. The most relevant quantities for representative 
x-strain values are listed on Table m 

Compressive x-strain of 4% slightly reduces Nq (see 
Table 0: and leads to a decrease in Tc to only 0.5 /iK. 
In contrast to this Nq, we see from Table 0 that Eccax 
increases with this compressive x-strain. Due to the pro¬ 
portionality between Emax and T *, the spin susceptibility 
starts to follow the logarithmic-T behaviour - the signa¬ 
ture of the VHS - at higher temperatures than the un¬ 
strained case. For example, we see that a compressive 
x-strain of 4% leads to T* of about 97 K (8.4 meV) [see 
Fig. I^a), yellow squares]. Due to the relatively high 
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Table II: Parameters related to the VHS for different 
strains on the zigzag direction [a:-axis in Figj^a)]. For 
T<T*, the susceptibility follows the logarithmic 
temperature dependence of Eq. 0. The PBEsol 
functional is used. 



^str — 49o 

Xatr = 0% 

2^str — +4% 

kmax (1/A) 

(0.175,0) 

(0.105,0) 

(0.033,0) 

Tmax (meV) 

27.8 

6.6 

0.1 

No (eV-i) 

0.0374 

0.0588 

0.1909 

T* (K) 

97 

17 

0.23 

Tc (K) 

5xl0~'^ 

4x10"® 

5x10"^ 

(a) 



(b) 




0 +1 +2 +3 +4 

x-strain (%) 

(C) 


N'o^eV) 


Eigure 4: (Color online) Effect of zigzag ridge [a:-axis in 
FigjPa)] strain on the VHS (a) Temperature dependence of 
the bare spin susceptibility y for different strains. The 
arrows indicate the temperature T* at which the 
susceptibility starts to deviate from the logarithmic 
behaviour, (b) Critical temperature as a function of 
strain, (c) Tc versus . The behaviour of Tc still follows 


the exponential law of Eq. (101 even if strain is applied. 


temperatures involved, the logarithmic-T behaviour in 
the spin susceptibility could in principle be observable 
experimentally, thus providing compelling evidence for 
the presence of the VHS. 

In contrast, tensile x-strain has the opposite effect 
on the band parameters: while Emax diminishes, Nq is 
greatly enhanced. Notably, the critical temperature ex¬ 
hibits an exponential dependence on tensile strain, as 
depicted in Eig. I^b). Eor a 4% strain, the critical 
temperature is about 0.05 K. Even though this corre¬ 
sponds to a 10^-fold increase in Tc with respect to the 
unstrained case, this magnetic state will hardly be seen 
experimentally due to the very low temperatures re¬ 
quired. We also observe that the logarithmic divergence 
becomes the dominant contribution at around T*~0.23 K 
(0.002 meV) [Fig. EJa), orange crosses], a much lower 
value compared to the unstained case resulting from the 


flattening of the valence band (and consequently dimin¬ 
ished Emax), caused by the applied stress. 

The physics behind the strain-dependence of the VHS 
is simple. The act of stretching will decrease the hoppings 
between phosphorus atomic orbitals, thus reducing the 
bandwidth. As the dispersion decreases, we expect that 
a and cJ in Eq. Q become smaller and hence the density 
of states Nq increases. The act of compressing will show 
the opposite trend. Because the band inflection is along 
kx, x-strain have a larger effect on a compared to y- 
strain, explaining our findings. 

Finally, we observe that the critical temperature still 
follows the exponential law of Eq. (fl^, even when strain 
is applied, as shown in Eig . [41 (c). For higher stress, 
Tc could deviate from Eq. ( |10[ ), since Emax diminishes 
(^10“"^ eV for a 4 % tensile strain) and therefore narrows 
the VHS divergence, limiting the increase of the critical 
temperature. 


IV. CONCLUSIONS 


We have used Wannier function-based interpolation 
techniques to investigate the VHS at the F point near 
the phosphorene Fermi energy with more than a billion 
fc-points. Thanks to this extreme resolution, we are able 
to present compelling numerical evidence for the presence 
of a VHS near the phosphorene Fermi energy. As a result 
of its close proximity to the valence band maximum, the 
VHS can be reached with a hole doping concentration 
on the order of 10^^ cm“^, easily achievable by chemical 
doping or ionic-liquid gating.l^ 

Furthermore, we have calculated an exact expression 
for the DOS near the VHS, and we have demonstrated 
that the spin susceptibility presents a logarithmic-T be¬ 
haviour, signature of the VHS, and consequent Liftshitz 
phase transition. 

We have also shown that the critical temperature can 
be increased up to 0.05 K by applying a modest strain to 
the phosphorene pliable waved structure. Although this 
ferromagnetic state would be very difficult to reach ex¬ 
perimentally, the logarithmic temperature behaviour of 
the spin susceptibility due to the presence of the VHS 
could be observed because it persists at higher temper¬ 
atures (T*^17 K for the unstrained case, and T*~97 K 
for a 4% tensile strain along the zigzag ridges). 

There are numerous experimental techniques able to 
detect the presence of VHSs. For example, the scan¬ 
ning tunnelling microscope (STM) measures the tun¬ 
neling differential conductance, which is proportional to 
the local DOS,^ and therefore represents an ideal tool 
to detect VHSs. This technique has been used to ob¬ 
serve VHSs in other 2D materials like twisted multilayer 
graphene,!^ or the cuprate superconductor Bi-2201.1^ 
Furthermore, angle-resolved photoemission spectroscopy 
(ARPES) can detect saddle points in the single-particle 
energy dis persion, a s employed for numerous cuprate 
compounds^iI^2lMIll and doped graphene.!^ Einally, the 




















Knight shiftP^in nuclear magnetic resonance experiments 
could provide evidence for the change in spin susceptibil¬ 
ity in proximity of the VHS. 
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Appendix A: Exact derivation of the density of 
states 


with 


k± = 


\ 



^max eJK 


(A3) 


Thus, if A < 0, k- is not purely real, and the lower 
bound is zero. Therefore, we have 



dk^ 




1.2 


-ikt-E 


1 

(-f) (.kl-kl){kl-k^_) 

-2i 1 

VPk-Jo y/il-x^)il-p^x^) 



-2i 


F 



-2i 


K{p). 


(A4) 


In this section, we present an analytical derivation of 
the DOS for phosphorene around the VHS. According 
to Eq. 0 . the dispersion relation of the valence band 
around the F point has the form 


where we have defined 


Ek = ^akl - ^f3k^ - ^a'k^ (AI) 


with a, (3, a' > 0. The valence band has its energy ex¬ 
treme at Emax = oP’jA.fd when (kx,ky)= (A:niax)0) and 
^max = \Jcilf3- Moreover, there is a VHS at energy 
A'vhs = 0 originating from states near kx = 0. 

By definition, the DOS per spin per unit area is 


NiE) = J 


dh dh 

^0^5{E-Ek) 


(A2) 


1 


dkx 


dk,. 


1 




dk. 


\dkyEk 

1 


-^iky - k^) 


“^2 _ _ z? 

2'^x 4 ^ 


where ky satisfies E = \akx — \[3k'^ — {ky^. 

The integration range is limited by the fact that the 
square root term has to be real. After some algebra, one 
can show that the integral range is fc S [max(0, fc-), fc+] 


and introduced F{(j), k) the incomplete elliptic integral of 
the first kind and K{k) the complete elliptic integral of 
the first kind. They are related here by F{^,k) = K(k). 

On the other hand, if A > 0 (fc_ > 0), the lower bound 
is k- and we will deal with 



-2 

''7W, 


-K{^Jl-p 2). 


In Eqs. (A4) and (A6|, we have used the integral formu¬ 
lae: 
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dx 


dx 


dz- 


V^(l - x‘^){l -p^x'^) 2 7q \/z{l - z){l -p^z) 


= F (arcsin(M),p) 


^ ^ =z!f( 

-\/(l - x‘^){l - p^x'^) 2 7„2 ip^z{l - z)(z -p~^) P y 


dz 


arcsin 


' 1-u^ 

1 — p~‘ 




(A7) 

(A8) 


For negative energies, E < 0, 

k- = i^\/l + \E\ /Amax = i\k-\ 


P = ^^l^ I 1 + Wl 


\E\ 


\E\ 

EmRTi 


(A9) 


= i \p\ (AlO) 


r 


and similarly, 


1 


N{E^0+) = —\ — 
ZTT^ V CtCt 


1 


In 


E„ 


E 


0 ( 1 ) 


(A16) 


After multiplication by the unit cell area axXciy, we ob¬ 
tain the result reported in Eq. Q. 


and therefore we will use the relation 


K{ik) = 


-.K 


ypTT V V fc2 -p 1 


fc2 


(All) 


in Eq. (A4). 


In conclusion, the density of states for E < 0 and E > 
0 are, respectively. 


N{E < 0) = 


1 




dkx 

1 


° ^J%kl-{kt-E 

K {i IpI) 

1 " ^ 


(A12) 


/3a' TT^i |fc_ 
1 


/ 3 «' ^^\f^-\^\pf + i^[\l^ + \p\ 


and 


N{E > 0) = 


dk 


fc- ^^kl-{kt-E 


(A13) 


2 1 
Pa' n’^k^ 


Vl-p- 


The final step is to analyze the asymptotic behaviour 
of the DOS. Since K will show a logarithmical divergence 
when 


K{k = l-rj) EP^^lii\p\ + 0(1), 


(AM) 


for if —)• 0, both the quantities ^ ^_|_|p |-2 and a/ 1 — p“2 
approach one. By using k± —)■ V^femax, |fc-| —t 


fcmaxY'2B!,L’ |/v| ^ -y |^| , 

we obtain the DOS at the VHS as: 


, H ^ 2.^/%-, and ^ 


N(E^0\ = ^ 

27r2 V aa 


In 


E„ 


\E\ 


+ 0 ( 1 ) 


(A15) 


Appendix B: Bare spin susceptibility at the VHS 
and the effect of doping 

Firstly, we derive the bare susceptibility when the 
Fermi energy is at ifvHS=0 (At=0). The spin suscepti¬ 
bility is given by 


X(T) = 




(27r)2 


1 

4T 


d^kE— cosh ^ ( -+ 


1 

'TT 

1 


rE„ 


' —oo 

/•A 


dEN/E) cosh 


^KjlT 


— 2 ( Ek 
TT 

E 

E 


=ivo r dxi.m^)cosi.-^x 
2 J-K/2T V 1^1 ) 

= — Nq dx 

Jq 

/A\ /■A/2T 

-f Vo In — 

° \2T 


cosh^ X 


dx cosh ^ X 


where we have considered only the contribution from the 
VHS and the logarithmic behaviour applies when |if | < A 
and A is another energy cutoff A<E^ax- 

If we then consider the limit A ^ T and use the fol¬ 
lowing formulae: 

= log Y — 7e = —O = —0.8188 (Bl) 


pOO 

/ 

/o cosh X 


( 7 e is the Euler-Mascheroni constant) 


dx cosh X = tanh — 

n 2il 


1 


(B2) 


we obtain the following expression for the bare suscepti¬ 
bility 

X(T) ^NoC + No ln( A) = No Hujd/T) (B3) 
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which is Eq. Q in the main text. 

Then, we discuss the effect of doping on the suscepti¬ 
bility. When we shift the Fermi energy from zero (Eyns) 
to fi, the VHS will change to —/r, and therefore we can 
replace the DOS by N(E) « Nq In for \E + fi\ < A. 

The susceptibility thus becomes 


A 


A 


.\E + f^\ 


cosh-^ f ^ 


(B4) 




dEluAcosh-^ 


We will now consider two regimes: T ^ fi and T <C /r. 
Let us start with the case By using the expan¬ 

sion 


cosh 


E — 
2T 


= cosh 




-1 

2TJ 


-anhl ALinh(A) 


(B5) 


and the approximation 

-2 f E — n 


cosh 


2T 


cosh-^ ( -] cosh-^ (A) 


(B6) 


the susceptibility can be written as 

rA 


X{T > /x) « AtVo y ^ dE In a cosh ^ ^ A^ (b7) 


X cosh ^ f—) 
\2TJ 

= No In (wd/T) cosh"^ 
« A^oln (ujd/T) 


(B8) 


and therefore, in this regime, the suscepti bilit y is the 
same as the undoped case [compare to Eq. (B3)]. 

On the other hand, when T<S^fj,, the function 
cosh“^('^^) decreases proportionally to exp(^^^^) for 
\E — ^\ ^ T. Therefore it is a good approximation to 
replace In cosh“^(^^|^) by In A cosh“^(^^|^) in the 


integrand in Eq. (B4). As a result, 

X{T < /i) « A A^^o y ^ dE In A cosh’ 


-2 f E- fj. 


1 A AT 

= — A^o / dEln-r—- cosh 

Jo \E\ _ 


AT 
+ cosh 
1 


2T 

-2 f E-fJ. 
2T 


+ 


-2 


E 




2T 

A A 

— iVo / dElr,- 
4T do M 

1 —2 I E ~\- fji 
+ cosh ( 


cosh 


-2 (E_-^ 
2T 


(B9) 


Then, using 


dE cosh-2 ( AA^ I = 2T 


2T 


tanh ( ± 


2T 




(BIO) 


we obtain 


X(T < ^J.) =^Noln^ 


tanh 


A 


2T 


' tanh 


A-/ 

2T 


sinh ^ 


= - Ah In — ■-T-;-7-- 

2 cosh cosh A# 


sA"o In — tanh — 
T 


' AToln 


2T 

A 


(Bll) 


The final result of Eq. pTj ) indicates that, in this 
regime, the bare spin susceptibility is independent of tem¬ 
perature, as seen in Fig. Sb). This behaviour originates 
from the infrared cutoff of the excitations due to the ther¬ 
mal distribution. Let us in fact consider the integrand in 
Eq. (B9|: 


In 


A 

\E\ 


cosh 


_2 /E-Zi' 


2T 


- cosh 2 


^ E + zx' 


A 


~ In —- cosh 2 
\E\ 


2T 

'E- 




(B12) 


2T 


Due to the chemical potential shift /x, the thermal dis¬ 
tribution function, cosh-2(AA^)~exp —^(A^)2 , is 
now centered at E=fi and not at E=0 like in the undoped 
case. Therefore, the integration around E=0, where the 
logarithmic function diverges, now makes essentially no 
contribution to the total integral, giving raise to the flat¬ 
tening of the susceptibility observed for T<Czx. 
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